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Recent marine heatwaves in the North Pacific 
warming pool can be attributed to rising 
atmospheric levels of greenhouse gases 


Armineh Barkhordarian® !™4) David Marcolino Nielsen’? & Johanna Baehr! 


Over the last decade, the northeast Pacific experienced marine heatwaves that caused 
devastating marine ecological impacts with socioeconomic implications. Here we use two 
different attribution methods and show that forcing by elevated greenhouse gases levels has 
virtually certainly caused the multi-year persistent 2019-2021 marine heatwave. There is less 
than 1% chance that the 2019-2021 event with ~3 years duration and 1.6°C intensity could 
have happened in the absence of greenhouse gases forcing. We further discover that the 
recent marine heatwaves are co-located with a systematically-forced outstanding warming 
pool, which we attribute to forcing by elevated greenhouse gases levels and the recent 
industrial aerosol-load decrease. The here-detected Pacific long-term warming pool is 
associated with a strengthening ridge of high-pressure system, which has recently emerged 
from the natural variability of climate system, indicating that they will provide favorable 
conditions over the northeast Pacific for even more severe marine heatwave events in the 


future. 
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(MHWs)—periods of anomalously high sea surface 

temperature (SST) at a particular location that lasts for at 
least 5 consecutive days—is projected to increase further in the 
twenty-first Century!~>. The probability of MHWs exceeding the 
pre-industrial 99th-%tile will increase under future global 
warming!. Attribution results by Laufkoetter et al.° further show 
that the occurrence probabilities of the duration and intensity of 
high-impact MHWs have increased 20-fold in the historical 
period (1982-2017) in comparison to the pre-industrial climate. 
MHWSs are extreme events with wide-ranging impacts on marine 
ecosystems, including geographical shifts of species’, and mor- 
tality of marine mammals and sea birds®. Beyond their impacts 
on marine ecosystems, ecological responses to MHWs can have 
socioeconomic implications’, such as loss of fisheries income, 
food availability, erosion of essential ecosystem services (e.g., 
carbon capture, water quality), and mass mortality of iconic 
species’. The aim of this study is to identify the fraction of the 
likelihood of a MHW event’s magnitude that is attributable to 
GHG forcing, which is important not only scientifically, but also 
for decision-making regarding ocean ecosystems, the economics 
of regional fisheries, and marine life®-!! as the GHGs emissions 
continue to rise. 

Over the last decade, the northeast Pacific ocean experienced a 
rapid resurgence of the Blob-like SST anomalies that caused 
devastating marine ecological impacts!!!*, such as very low- 
ocean productivity!’; dramatic mortality events in seabird 
species!4; and major outbreaks of harmful algal blooms that 
produce extremely dangerous toxins!®. Furthermore, they con- 
tributed to severe drought conditions across the US West Coast!®. 
The 2014-2015 MHW!7-!°, the so-called first warm “blob”?%7!, 
which appeared off the coast of Alaska, originated in winter and 
was the result of a resilient atmospheric ridge in the Northeast 
Pacific, which weakened the climatological Aleutian Low and 
related surface winds*®. The prolonged nature!® of the 
201-2015 MHW is suggested | to be linked to the dynamics and 
coupling between the two dominant modes of winter 
(January-March) SST variability in the North Pacific, namely the 
Pacific Decadal Oscillation (PDO2*), and the North Pacific Gyre 
Oscillation (NPGO?3). A stronger NPGO-PDO coupling is pre- 
dicted under anthropogenic forcing!’. The second Blob (2.0) 
peaked in summer 2019 and resulted from a weakened North 
Pacific High, which reduced the strength of the surface winds, 
resulting in reduced evaporative cooling in the Northeast 
Pacific?+?°, The study by Holbrook et al.!? identifies multiple 
drivers that can enhance the occurrence of MHWs over the North 
Pacific, including El Nifio and the positive PDO phases. 

Unlike previous studies which have focused on linking the SST 
patterns in the North Pacific to changes in the oceanic circulation 
and the extratropical/tropical teleconnections”!?-17:1820,24,26, we 
here perform two different statistical attribution methodologies in 
order to identify the human fingerprint in Northeast Pacific SST 
changes both on multidecadal timescale (changes of mean SST) 
and on extreme SST events on daily timescale (Marine Heat- 
waves). Evidence that anthropogenic forcing has altered the base 
state (long-term changes of mean SST) over the northeast Pacific, 
which is characterized by strong low-frequency SST fluctuations, 
would increase confidence in the attribution of MHWs7’, since 
rising mean SST is the dominant driver of increasing MHW 
frequency and intensity, outweighing changes due to temperature 
variability! 

In this study, we provide a quantitative assessment of whether 
GHG forcing, the main component of anthropogenic forcings, 
was necessary for the North Pacific high-impact MHWs (the 
Blob-like SST anomalies) to occur, and whether it is a sufficient 
cause for such events to continue to repeatedly occur in the 


O n the global scale, the frequency of marine heatwaves 


future. With these purposes, we use two high-resolution observed 
SST datasets, along with harnessing two initial-condition large 
ensembles of coupled general circulation models (CESM1-LE?®° 
with 35 members, and MPI-GE*? with 100 members). These large 
ensembles can provide better estimates of an individual model’s 
internal variability and response to external forcing?!3*, and 
facilitate the explicit consideration of stochastic uncertainty in 
attribution results*>. We also use multiple single-forcing experi- 
ments from the Detection and Attribution Model Intercompari- 
sion Project (DAMIP*4) component of Coupled Model 
Intercomparison Project phase 6 (CMIP6°°). 

Two statistical methodologies are used: (1) with the first 
method, we analyze the observed long-term spatiotemporal 
changes of SST to (a) detect the presence of a signal beyond 
changes solely due to natural (internal) variability, and to (b) 
attribute the detected changes in long-term SST to external cli- 
mate drivers. The climate over North Pacific is potentially 
influenced by two external climate drivers: well-mixed green- 
house gases (GHGs) and anthropogenic aerosols (AER), which 
have opposing effects on the radiation energy balance*®. Our 
attribution analysis is based on assessing the amplitude of the 
response of SST to each external forcing from the observations via 
the estimation of scaling factors>”>°. (2) The second statistical 
method is extreme event attribution?”*°, which determines how 
anthropogenic forcings have changed the likelihood of occurrence 
of a particular event. Following Hannart et al.4°4!, we present 
extreme event attribution in terms of necessary and sufficient 
causation. Our results, based on the two different attribution 
methods, provide a complete picture of anthropogenic influence 
on SST extremes over the North Pacific. 

In this study, we show that forcing by elevated well-mixed 
GHG levels, has virtually certainly caused the multiyear persistent 
2019-2021 marine heatwave, in a necessary causation sense. 
There is less than 1% chance that the 2019-2021 event with ~3 
years duration and 1.6 °C intensity could have happened in the 
absence of GHG forcing. We further discover an outstanding 
warming pool over the Northeast Pacific, co-located with the 
multi-year persistent MHW events. The warming pool is marked 
by concurrent and pronounced increase in annual mean, and 
variance of SSTs, decrease in cold-season low-cloud’s cooling 
effect, and strengthening cold-season atmospheric ridge of high- 
pressure system. There is less than 5% chance that natural 
(internal) variability is responsible for the observed cold-season 
strengthening high-pressure system. Our long-term SST detection 
and attribution results further reveal that forcing by elevated 
GHG levels and the recent industrial aerosol-load decrease are the 
key causes for the configuration of the here-detected long-term 
warming pool (P< 0.05). 


Results 

Marine heatwaves characteristics. We detect 40 MHW events 
over the Northeast Pacific from January 1982 to March 2022 in 
NOAA OISSTv2* daily SST data (Fig. 1). The detected MHWs 
over the twenty-first century (2000-03/2022) are 4.5-fold more 
frequent, ninefold longer-lasting, and threefold more intense in 
comparison with those occurred in the previous decades (Fig. Ic). 
During the three major MHWs over the northeast Pacific in 
2014-2015 (Fig. 2b), 2016 (Fig. 2c), and 2019-2021 (Fig. 2a) the 
maximum SST anomalies (above the climatology in the peak date 
of the event) reached 5.6 °C, 5°C, and 5.8 °C, respectively. The 
2014-2015 MHW, which received major societal concerns, lasts 
for 600 days with 1.4°C intensity (average SST anomaly during 
the event). The 2019-2021 MHW, which lasts for almost 3 years 
(1000 days) with 1.6 °C intensity, is the most severe MHW both 
in terms of intensity and duration that has ever been detected 
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Fig. 1 External forcing had a significant contribution to the observed marine heatwaves characteristics. a Intensity, b duration of the observed 40 marine 
heatwave events (MHWs) detected between January 1982 to March 2022 over the northeast Pacific (160°W-100°E, 10°N-60°N). Each black (red) bar 
represents a MHW event and the gray (blue) whiskers show the 98%-tile distribution of MHW characteristics in an undisturbed Pre-industrial climate. 
Detection is claimed in cases where the whiskers do not include the zero line (Ho: Obs. + Pog, # 0). Events that are statistically distinguishable (P < 0.02) 
from the distribution of events in Pre-industrial climate are denoted by red bars with blue whiskers. Up to 60% of the events detected over the last decade 
(2010-2021) are either more intense and/or longer-lasting than could solely be attributed to climate variability in the absence of external climate drivers. 
The two major MHWs with high impact occurred in 2014-2015, and 2019-2021 are denoted with stars. e Changes in MHWs characteristics (duration, 
frequency, and intensity) over 2000-03/2022 in comparison with those detected over 1982-1999. 


Evolution of the MHW from 2019 to 2021 b Evolution of the MHW from 2014 to 2015 
O°, 
July 2019 November 2019 July 2020 . °c Sarena 2085 September 2014 June AOtS a 2015 C > 
a 15 @ os & 
as 5 oo 3 
=08 hed EF io 
Bi 8 3s B 
-25 
September 2020 November 2020 August 2021 °C Evolution of the MHW from 05/2016 to 11/2016 
* ‘ : , [: 8 C August 2016 September 2016 °C 
1 ae 2 iE, . 
Hos § - x a (1%. 3 
H-05 asl 1 € 
I: a i 05 6 
15 no if q . 0 (> 
‘ 2 05 © 
-25 J -1 oa 
15 
3s 
Annual mean SST linear trend SST variance linear trend Cold-season 500Gph linear 
(OISST; 1996-2021) € — (OISST; 1996-2021) f trend (ERAS; 1996-2021) 
» ; 0.4 7 ee Tica = Pow 30 
8 0.3 8 r 8 4 es 
0.2 0 15 @ 
z o1 8 z 8 z wv 0 
= 0 8 g 3) 8 5 8 
® 2 oe ® 
= 0:1 aa 2 an ra -10 z 
8 -0.2 9 g i) 8 15 € 
~20 
-0.3 -25 
_-e " -0.4 : emai i 
“90°E 150°E 150°W 90°W 90°E 150°E 150°W 90°W 


Fig. 2 Evolution of the MHWs, and the long-term trends of the background state. a SST anomaly patterns (above the 1983-2012 climatology) during the 
evolution of the MHW from 2019 to 2021, b the MHW from 2014 to 2015, and ¢ the MHW from 05/2016 to 11/2016. Observed pattern of trends over 
1996-2021 in (d) annual mean SST (OISST; Units: °C decade—!), e SST annual variance (OISST; units: °C2 decade), f cold-season (November-to- 
February) geopotential height at 500-hPa (ERAS reanalysis; units: m decade—!). 
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Table 1 Attribution of MHWs duration and intensity to GHG forcing. 


Date of event 
2014-2015 

2016 
2019-2021 


Threshold of intensity 
1.4°C 
11°C 
1.6 °C 


Threshold of duration 


600 days >10& 
163 days 2.4 
1000 days >10& 


Date of event 
2014-2015 


2016 
2019-2021 


PR of intensity 


R of duration 


PN of intensity 
1.0 [0.98-1.0] 
1.0 [0.97-1.0] 
1.0 [0.98-1.0] 


PS of intensity 
0.00 
0.05 
0.00 


PS of duration 
0.02 


0.08 
0.02 


PN of duration 
1.0 [0.98-1.0] 


0.58 [0.57-0.78] 
1.0 [0.98-1.0] 


We present the event-attribution results as the probability of necessary (PN) and sufficient (PS) causation as well as the probability ratio (PR) for three high-impact MHWs detected over the North 
Pacific, 2014-2015, 2016, and 2019-2021. 


since the year 1982. This raises the question of whether their 
occurrence, intensity, and duration have been influenced by 
external climate drivers. 

Our results reveal that up to 60% of the MHW events detected 
over the last decade (2010-03/2022) are either more intense 
(Fig. 1a) and/or longer-lasting (Fig. 1b) than could solely be 
attributed to climate variability in the absence of external climate 
drivers. We arrive at this conclusion by comparing the detected 
MHWs with events occurring solely in response to the natural 
variability of the climate system. To obtain these estimates, we use 
CMIP6 pre-industrial control simulations of daily SST, which 
provide pseudorealizations of MHWs characteristics in the 
absence of external climate drivers#?*4. In this manner, we can 
identify events that are significantly distinguishable from the 
distribution of events detected in Pre-industrial climate (Ho: 
Obs. + Pogoy x 0). 

Here, we conclude, therefore, that external forcing has a 
significant contribution (P < 0.02) to the intensity and duration of 
the detected MHWs, including the three major 2014-2015, 2016, 
and 2019-2021 MHWs. However, it remains unclear whether the 
well-mixed GHG forcing is a sufficient causation for the 
occurrence of these MHWs and to what extend GHG forcing 
has changed the likelihood of the observed events. In the 
following, we use “Extreme Event Attribution?”**” framework to 
identify the fraction of the likelihood of these events duration and 
intensity that is attributable to GHG forcing. 


Marine heatwaves attribution (extreme event attribution). We 
use extreme event-attribution technique*’, and assess the influ- 
ence of GHG forcing on the duration and intensity of the 
observed MHWs. In extreme events attribution approach?”4>*6, 
Y is defined as an observed extreme situation based on excee- 
dance over a threshold u of a relevant climate index Z. We 
evaluate the extent to which an external climate forcing f, for 
instance greenhouse gas (GHG) forcing, has changed the prob- 
ability of occurrence of the event Y. The variable X; is introduced 
to indicate whether or not the forcing fis present. The probability 
pl=P(Y=1|X;=1) of the event occurring in the real world, 
with f present, is referred to as factual, while pp = P(Y = 1|X;= 0) 
is referred to as counterfactual (with forcing f being absent)*>*°. 
Following Hannart et al.404!, we present event attribution in 
terms of necessary and sufficient causation. According to their 
definitions, requiring the presence of a particular forcing scenario 
(f) for an event to occur would be necessary causation. In con- 
trast, if the particular forcing scenario (f) always produces the 
event in question, this would be sufficient causation*?*!. In this 
study, the event is MHW’s intensity and duration (Y) and that 
particular forcing scenario (f) is GHG forcing. 

For each detected MHW, we estimate the probabilities that a 
MHW has occurred that equals or exceeds the duration and 


intensity of the observed MHW in ALL forcing (actual world) 
and fixGHG forcing (counterfactual world) scenarios. That is, we 


calculate the probability, EN GHIG > of the threshold (duration of 


the observed MHW) being exceeded without GHG forcing, and 
the probability, P#“"*°", of exceeding the threshold with GHG 


forcing. Similarly, we calculate the probability, Pe , of the 


threshold (intensity of the observed MHW) being exceeded 


without GHG forcing, and the probability, Pi7y"”, of exceeding 
the threshold with GHG forcing. Choosing the 1982-2021 time 
period, the MHW characteristics (duration and intensity) from 
each year and each 20 realizations of CESM1 model are pooled 
together (780 years total) and probabilities are estimated. 

The estimated probabilities are used to calculate three event- 
attribution metrics*°4” (see “Methods”). The probability ratio 
(PR), which describes how many times as likely the event 
occurrence is with ALL forcing than ALL minus GHG forcing 
(fixGHG). The probability of necessary causation (PN) describes 
the probability that GHG forcing is a necessary cause of the 
particular event; that is, that GHG forcing is required for the 
event’s occurrence, and finally the probability of sufficient 
causation (PS) describes the probability that GHG forcing is 
sufficient for the event, such that a scenario with GHG forcing 
will see the occurrence of this event every time. A summary of the 
PS, PN, and PR values for the duration and intensity of the three 
selected MHWSs is presented in Table 1. 

The probability ratio values (PR, Eq. (1), Eq. (2)), which 
quantifies the additional probability of an event’s intensity 
(duration) due to GHG forcings, increases for the more extreme 
MHWs (Table 1). A PR value of 10° for both 2014-2015 and 
2019-2021 MHWs implies that the occurrence of such events is 
10° times more likely under the influence of GHG forcing. The 
probability of necessary causality (PN, Eq. (1), Eq. (2)), which 
describes the probability that GHG forcing is a necessary cause of 
the detected MHWs, increases with increasing the severity of 
MHWs in terms of both duration and intensity. For instance, a 
PN value of approximately 0.99 ([0.98-1.0]) for the 
2019-2021 MHW means that 99% of the probability of such an 
event is due to GHG forcing. In other words, there is a 99% 
chance that GHG forcing is required for this event to occur. The 
probability of sufficient causation (PS, Eq. (1), Eq. (2)), which 
describes the probability that the inclusion of GHG forcing is 
sufficient for the event’s occurrence, is small (<0.01) for the more 
extreme events (Table 1), as such events are rare even with ALL- 
forcing scenarios. 

In summary, event-attribution results provide evidence that the 
occurrence of a multiyear persistent MHW, such as that occurred 
in 2014-2015 (with 600 days duration and 1.4 °C intensity) and in 
2019-2021 (with 1000 days duration and 1.6°C intensity), are 
entirely attributable to the combination of anthropogenic and 
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Fig. 3 Detection of externally forced trends in large-scale circulation patterns. a Observed trend pattern in cold-season (November-to-February) over 
1996-2021 of 500-hPa geopotential height (500 Gph; units: m decade~'), and b of mean sea-level pressure (SLP; units: Pa decade—'). Regions where 
systematically and externally forced changes are detectable at 95% confidence level in (¢) 500 Gph, and (d) SLP, in comparison with 280 
pseudorealizations of unforced trends in 500 Gph and SLP due to pre-industrial variability derived from CMIP6 control runs. e Observed pattern of trends 
over 1995-2018 in total cloud-cover fraction (EUMETSAT; units: % decade—'). f Changes in daily mean SST and seasons based on NOAA OISST over the 
Northeast Pacific warming pool (indicated with a black box on Fig. 4c). Expansion of summers by more than one month and shrinking of winters over the 


warming pool. 


natural forcings (ALL forcing). The occurrence of the MHWs was 
extremely unlikely in the absence of GHG forcing (with less than 
1% occurrence probability under no-GHG effect), and GHG 
forcing is necessary for the occurrence of these events (PN = 1.0 
[0.98-1.0]; Table 1). Those MHWs are extreme in the current 
climate, so the inclusion of GHG forcing is a necessary, but not 
just merely sufficient, causation (PS <0.01; Table 1). In other 
words, a MHW of this magnitude requires GHG forcing to occur 
(with >99% probability), but the inclusion of GHG forcing alone 
is not enough to guarantee the event’s occurrence. 


Discovering an outstanding long-term warming pool. We dis- 
cover an outstanding warming pool over the Gulf of Alaska in the 
high-resolution NOAA OISST*4? dataset, co-located with the 
MHWs induced SST anomalies. The region is marked by a con- 
current and pronounced increase in annual mean (0.4 °C decade“!; 
Fig. 2d), and annual variance of SSTs (0.6 °C? decade~!; Fig. 2e) 
over 1996-2021. Hereafter, we refer to this region as the Pacific 
long-term warming pool. 

The co-location of the major MHWs over the Northeast Pacific 
(Fig. 2a-c) with the here-detected long-term warming pool points 
towards possible positive feedback, and suggests that the 
prominent MHWs occur where mean warming (Fig. 2d), and 
higher variance overlap (Fig. 2e). We also detect significant 
changes in daily mean SST and seasons over the warming pool. 
SSTs are now (2000-2021) higher on average for every day of the 
year than they were in the late twentieth Century (1982-1999). 
Over the last two decades, summers are on average 1.5 °C warmer 
and 37 days longer. SSTs that marked the start of summer (June 
1) now come 11 days earlier, and SSTs that marked the end of 
summer (September 1) now come around 27 days later. Winters 
are on average 0.5 °C warmer and 11 days shorter (Fig. 3f). 


Large-scale atmospheric circulations play an important role in 
the occurrence of extremes*®. In the cold season (November-to- 
February), the observed trends in geopotential height at 500-hPa 
(500 Gph) according to ERA5*° reanalysis show a tendency 
towards a ridge with a magnitude on the order of ~+30m 
decade! increase, located over the warming pool area (Fig. 3a), 
resembling the ridiculously resilient ridge®°>!, which we show that 
is getting intensified over the 1996-2021 time period. The cold- 
season enormous strengthening ridge is associated with a 
~+250 Pa decade! increase in mean sea-level pressure (SLP, 
Fig. 3b), resulting in an amplification of the background state of the 
MHWs. The pronounced increase in 500Gph and SLP, which 
enforce heat-trapping systems, are also associated with a decreasing 
trend in cloud cover starting in 1995/1996. The EUMETSAT** 
satellite data shows a 5% decade! decreasing trend in cold-season 
cloud cover during 1995-2018 (Fig. 3e). Low-cloud cover reduction 
is the major contribution to the observed decline in total cloud 
fraction (not shown), resulting in decreases of winter-time low- 
cloud’s cooling effect. 

The detected long-term pacific warming pool could have major 
ecological consequences. For instance, marine species shift their 
distributions in response to warming trends”!!, and the size of 
various fishes decline due to increasing water temperature and 
reduced oxygen levels°*>4, causing changes to the entire 
ecosystem. 


Detection of systematic changes in the Pacific warming pool. 
However, it remains unclear whether the here-detected Pacific 
long-term warming pool, and the associated trends in large-scale 
circulation patterns (SLP and 500 Gph), are a simple manifestation 
of internal climate variability or are externally and systematically 
forced. To address this question, we begin with comparing the 
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observed trends with estimates corresponding to the natural 
(internal) variability of the climate system, which ENSO, PDO, 
NPGO etc. are part of. To obtain estimates of the natural (internal) 
variability we use three data sources (see “Methods”). We test the 
null hypothesis that the observed trends in SST, SLP, and 500Gph 
over 1996-2021 are within the 2.5-97.5%-tile distribution of 
unforced trends of those variables (as derived either from the pre- 
industrial climate or historical variability) or naturally forced 
trends (as derived from the 850-1850 millennium simulation)3*°>. 

There is less than 5% chance that internal variability is 
responsible for the observed cold-season strengthening high- 
pressure system in the region. In spite of the large interannual 
variability over the north Pacific region, our detection analyses 
(displayed in Fig. 3c, d) show that the pronounced increases in 
500 Gph (+30m decade~!) and SLP (+250hPa decade™!), 
respectively, over 1996-2021 are systematically forced and larger 
than could be due to natural (internal) variability alone (with 
<0.05 risk of error). Such changes serve to systematically enhance 
atmospheric stability, supporting the configuration of the 
warming pool. 

Detection analyses on the SST trends (displayed in Fig. 3c-e) 
reveal that the SST increase over the warming pool is as well 
stronger than could be due to natural (internal) variability alone, 
and systematically forced changes of SST are detectable at the 
95% confidence level in summer and autumn (JJASON; 
June-November). This result is robust across comparisons of 
observed SST trends with (1) unforced trends derived from 
CMIP6*° pre-industrial control simulations, which provide up to 
280 pseudorealizations of how SST might have changed in the 
absence of external climate drivers (Fig. 4b), (2) trends obtained 
from time-evolving historical variability over 1996-2021, derived 
from the MPI-ESM GE?9 100-member ensemble (Fig. 4c), and (3) 
naturally forced trends derived from CCSM4°° Paleo simulations 
over 850-1850 (Fig. 4d). 

On the one hand, in winter (JDF), and in spring (MAM) a 
substantial portion of SST trends over the North Pacific can be 
explained by the natural variability alone, over 80% and 75% of 
the grids, respectively (Supplementary Fig. 1). On the other hand, 
annually, in summer (JJA), and in autumn (SON), the 
configuration of the Pacific warming pool emerges from the 
natural variability and is found to be externally and systematically 
forced (P < 0.05). 

Having demonstrated that the here-detected warming pool 
significantly deviates from natural variability, that is, externally 
forced changes are detectable over the region (P<0.05), we 
proceed to an attribution step by checking whether the detected 
changes are consistent with what climate models describe as the 
expected response to anthropogenic forcing. The following 
attribution analysis focuses on JASON months, where the most 
significant and pronounced warming trends in SST are detectable 
over the last decades (Supplementary Fig. 1). 

As a first step to determining if the observed SST trends 
constitute the system’s forced response to external climate drivers, 
we project the observed SST changes on the model simulated 
forced response, ALL signal (anthropogenic + natural) by using 
univariate total least square (TLS) regression analysis (Eq. (3) in 
“Methods”). The forced response (ALL signal) is derived from 
two ensembles of model simulations: the CESM1 35-member 
Large Ensemble (CESM-LE?°) and the MPI-ESM 100-member 
Grand Ensemble (MPI-GE?’). 

The 1-year moving scaling factors (a;) during 1983 to 2021 
time period demonstrates that the 25-year SST trends in the 
warming pool emerge from natural variability in 2018 and later 
on (Fig. 4e). We reach this conclusion, as the 95%-tile range of 
the internal variability-generated uncertainty of scaling factors for 
the row (dark gray), and double (light gray) the model variance 


(assessed from fits of the regression model (Eq. (3)) to control run 
segments) does not include the zero line but is consistent with 
unity (a, # 0 () a; = 1) for 25-year trends ending in 2018 and later 
on (1993-2018, 1994-2019, 1995-2020, and 1996-2021; Fig. 4e). 
This indicates the emergence of a detectable ALL signal 
(anthropogenic forcing being dominant) in SST changes in the 
twenty-first century (with <5% risk of error) over the warming 
pool. The similarity of the forced response in CESM-LE (Fig. 4e) 
and MPI-GE (Supplementary Fig. 2) suggests that the smaller 
ensemble size of the CESM-LE does not affect the forced signal. 

It is interesting to note that, the range of internal variability- 
generated uncertainty in the scaling factors (the gray shaded area) 
decreases when more recent years are included in the analysis 
(Fig. 4e). This could imply that the signal of SST increase is 
becoming dominated by the external forcing, which makes it 
easier to separate the signal’s large contribution from the internal 
climate variability (noise). In other words, the region might 
become more responsive to anthropogenic forcing as a result of, 
for instance, decreasing mixed layer depth that could cause a 
stronger SST response for the same heat flux”. 

The climate over North Pacific is potentially influenced by two 
external climate drivers: well-mixed greenhouse gases (GHGs) 
and anthropogenic aerosols (AER), which have evolving con- 
tribution in driving SST over the Gulf of Alaska. In the following, 
we focus on 1996-2021 time period and assess the response of 
SST to the GHGs-only forcing and the AER-only forcing in 
isolation. 


Evolution of forced SST trends. Here we are isolating the 
evolving contributions of industrial aerosols, AERindust (Note that 
biomass burning aerosols are not considered), and greenhouse 
gases (GHGs) in driving annual SST trends during the 1920-2021 
time period (Fig. 5a). We use two ensembles that are identical to 
the CESM-LE?®3?, but each ensemble excludes the time evolution 
of one forcing agent: greenhouse gases (LE-fixGHG, include 20 
members) and anthropogenic aerosols (LE-fixAERindust include 
20 members). We define the response to the withheld forcing as 
ie, ALL. — fixGHG,,,22. The time series of GHG-induced SST 
anomalies shows a steep positive trend, with an increase in 
magnitude over 1920 to 2025 time period (Fig. 5a). The SST trend 
pattern driven by AER;ndusi-forcing only shows that AERindust 
cools North Pacific over the 1980’s (Fig. 5a), and offsets GHG- 
induced radiative warming (Fig. 5b), with the most pronounced 
SST declines over the Northeast Pacific during 1956-1980 
(Fig. 5c). Noteworthy, aerosols contribute negatively to the 
radiative forcing response (cooling effect), and there is indeed an 
observed increase in aerosols in this period*®. However, in the 
most recent period, the cooling due to AER-forcing is replaced by 
warming (Fig. 5a, g). The reversal of the cooling trends due to 
AER-forcing to warming after around 1995/1996 is very likely a 
result of concurrent aerosols optical depth (AOD) declines over 
North America and Europe®®°?. Consequently, the ALL-forcing 
trend pattern over 1996-2021 is characterized by amplified 
warming, when both GHG-induced radiative warming and 
AER ndust-induced local warming have contributed, individually, 
to the SST increase (Fig. 5e). 

Interestingly, the AERindust-forcing only simulation over 
1996-2021 shows a horseshoe-like pattern of SST warming that 
is reminiscent of the PDO (Fig. 5g). Despite the PDO being 
commonly recognized as an internally generated mode of the 
climate system, recent studies have suggested that PDO phase 
transition is also significantly affected by external climate 
drivers®°-®?. Regional trends in anthropogenic aerosols over 
North Pacific, with aerosols-cloud interaction being the dominant 
aerosols forcing®’, can influence the PDO through modulation of 
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Fig. 4 Detection of systematically forced trends in observed sea surface temperature (SST) record. a Observed trend pattern of SST over 1996-2021 in 
JJASON (June-to-November) based on OISSTv2. Regions where systematically forced changes of SST are detectable at 95% significant level (b) in 
comparison with 280 pseudorealizations of unforced trends of SST due to pre-industrial variability derived from CMIP6 control runs, (¢) in comparison 
with trends due to historical variability over 1996-2021 estimated from the 100 realizations of MPI-GE and, (d) in comparison with 39 pseudorealizations 
of naturally (internal + external) forced trends derived from CCSM4 850-1850 Paleo simulation (P < 0.05). e One-year moving scaling factors of observed 
25-year SST trends onto the ALL signal (anthropogenic + natural) derived from the ensemble mean of CESM1-LE 35 members over 1983 to 2021 in 
JJASON over the warming pool (black box in ¢). The gray shaded area displays the 95%-tile range of the internal variability-generated uncertainty of 
scaling factors (a) assessed from fits of the regression model to 280 control run segments for the raw (dark gray) and double (light gray) the model 
variance. Detection of ALL signal is claimed when the gray shaded area does not include the zero line but is consistent with unity (a; # 0 () aj = 1, with <5% 


risk of error). 


the Aleutian low®?-°®, The CESMI1-LE model simulates a 
strengthening of the Aleutian Low in DJF (decrease in mean 
sea-level pressure in the north Pacific index (NPI) region defined 
by the 160-220°E, 30-65°N latitude-longitude domain) in 
response to AER-forcing over 1996-2021 (Fig. 5h). The winter- 
time strengthening of the Aleutian Low (—50Pa decade—! 
decreasing trend in SLP), induce SST trends in the northeast 
Pacific that resemble the PDO (Fig. 5g). This makes it difficult to 
distinguish PDO-related anomalies from climate change signals. 
Thus, it is plausible that the anthropogenic aerosols have 
influenced both on North Pacific SSTs and the PDO index itself. 


Univariate trend detection and bivariate trend attribution. The 
univariate detection analysis, which is based on estimating the 
amplitude of the response of SST to each external forcing from 
the observations via the estimation of scaling factors*”>*°7 (see 
“Methods”), reveals that irrespective of the models used, the 
elevated GHG concentration has a robust detectable influence on 
the observed SST increase over the region (P< 0.05, Fig. 6a). As 
the internal variability-generated uncertainty of scaling factors (a; 
in Eq. (3)) does not include the zero line for the GHG signal 
derived either from CESMI1-LE (ensemble mean _ of 
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20 simulations, green bar, Fig. 6a) or CMIP6-DAMIP models 
(ensemble mean of 34 simulations conducted by 5 DAMIP 
models, cyan bar, Fig. 6a). Neither the scaling factors of the AER 
signal derived from CESM1-LE (red bar) nor from CMIP6- 
DAMIP models (purple bar) include the zero line, suggesting that 
the AER-induced local warming over the Gulf of Alaska also 
contributes significantly to the observed increase in SST over the 
region in JJASON. 

Furthermore, in order to separate the driver’s contributions to 
the response, a combined influence of GHG, and AER should be 
considered. In this case, a bivariate (two-dimensional) attribution 
analysis is required, where the observed trend of SST is projected 
onto two hypothetical signals (GHG and AER) simultaneously”. 
The two-dimensional (bivariate) uncertainty contour for the GHG 
and AER is shown with an ellipse (Fig. 6c). The ellipse containing 
90% of the estimated joint distribution of scaling factors for the two 
signals (GHG and AER) excludes the origin (0,0), indicating that 
the effects of GHG and AER signals are detectable simultaneously. 
In addition, the scaling factors are consistent with unit amplitude 
since the point (1,1) lies within the ellipse (Fig. 6c). 

In summary, we detect a signal from both GHG and AER in 
the observed long-term trend of SST in JJASON, and demonstrate 
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Fig. 5 Evolving roles of the greenhouse gases forcing and industrial aerosols in driving annual SST over the North Pacific. a Area mean annual SST 
anomalies over the Gulf of Alaska (black dotted box in e), derived from the ensemble mean of 20 realizations of CESM1-LE model, ALLem minus fixGHGen 
in red, and ALLem minus fixAERem in blue. Response of annual SST to the ALL forcing (anthropogenic + natural), greenhouse gases (GHG) forcing and 
industrial aerosols (AER) forcing (b-d) over 1956-1980, and (e-g) over 1996-2021. h Response of mean sea-level pressure (SLP) to AER-forcing in DJF 
(December-February). The Aleutian low region defined by 160-220°E, 30-65°N latitude-longitude is denoted by a red box in h. 
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Fig. 6 Detection and attribution of sea surface temperature trends over the warming pool. a Univariate attribution: Scaling factors (a) of observed SST 
changes over the warming pool (black box in b), and its 95-%tile uncertainty range, assessed from fits of the regression model to control run segments, in 
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conducted by 5 models; cyan and purple bars, respectively). b Observed trend pattern of SST over 1996-2021 time period. ¢ Bivariate attribution: The ellipse 
displays the 90% of the estimated joint distribution of scaling factors for the GHG & AER signals when observed data are regressed onto two signals 
simultaneously during 1996-2021. The one-dimensional uncertainty intervals for the univariate and bivariate analysis for two signals are shown as red and black 
whiskers, respectively. Bivariate attribution is claimed in cases where the ellipse excludes the origin (0,0) but the point (1,1) lies inside the ellipse. 


8 COMMUNICATIONS EARTH & ENVIRONMENT | (2022)3:131 | https://doi.org/10.1038/s43247-022-00461-2 | www.nature.com/commsenv 


COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-022-00461-2 


ARTICLE 


that AER-forcing can act in favor of GHG-induced warming, 
amplifying this warming. Therefore, we conclude that forcing by 
elevated GHG levels and the recent industrial aerosol-load 
decrease are identified as key causes for the here-detected long- 
term warming pool (P< 0.05). 


Conclusions. Over the last decade, the North Pacific experienced 
strong marine heatwaves (MHWs) that produced devastating 
marine ecological as well as socioeconomic impacts, and received 
major societal concerns. We use the extreme event-attribution 
technique based on causal counterfactual theory*®4!47, and 
provide a quantitative assessment of whether GHG forcing, the 
main component of anthropogenic forcings, was necessary for the 
North Pacific high-impact MHWs (the Blob-like SST anomalies) 
to occur, and whether it is a sufficient cause for such events to 
continue to repeatedly occur in the future. 

In this study, we show that forcing by elevated well-mixed 
GHG levels, has virtually certainly caused the multiyear persistent 
2019-2021 marine heatwave. In other words, there is less than 1% 
chance that the 2019-2021 event with ~3 years duration and 
1.6°C intensity could have happened in the absence of GHG 
forcing. Our extreme event-attribution analysis further reveals 
that the GHG forcing is a necessary, but not sufficient, causation 
for the multiyear persistent MHW events in the current climate, 
such as that happened in 2014-2015, and 2019-2021. That is, a 
MHW of this magnitude requires GHG forcing to occur (with 
99% probability), but the inclusion of GHG forcing alone is not 
enough to guarantee the occurrence of the events. However, given 
that the occurrence of the 2019-2021 (2014-2015) MHW was 
extremely unlikely in the absence of GHG forcing (with <1% 
occurrence probability under no-GHG forcings), combined with 
increasing trends in MHWs duration and intensity, it is likely that 
future MHW events will be attributable to GHG forcing, and that 
the inclusion of GHG forcings will become a sufficient cause for 
events of the magnitude of the 2019-2021 (2014-2015) 
record event. 

We further discover a systematically-forced outstanding warm- 
ing pool in the high-resolution NOAA OISST* dataset, co- 
located with the MHWs maximum SST anomalies. The region is 
marked by concurrent and pronounced increase in annual mean 
(0.4°C decade), and variance of SSTs (0.6°C? decade), 
decreases in winter-time low-cloud’s cooling effect, and increases 
in atmospheric stability with a strengthening ridge. The region is 
further identified with 0.5°C warmer and shorter winters with 
37 days longer summers. Consequently, greater exposure to heat, 
and lack of usual winter-time cooling leads to 4.5-fold more 
frequent, ninefold longer-lasting, and threefold more intense 
MHWs in the twenty-first century, in comparison with those 
detected in the twentieth century. We further show that up to 
60% of the MHW events detected over the last decade 
(2010-2021) are either more intense and/or longer-lasting than 
could solely be attributed to climate variability in the absence of 
external climate drivers (P< 0.05). We refer to this region as the 
Pacific long-term warming pool. 

The univariate detection analysis indicate that anthropogenic 
signal has recently emerged from the natural variability of SST 
over the warming pool in JJASON (June-to-November). The 
pronounced increases in 500 Gph (+30m decade!) and SLP 
(+250 Pa decade—!) in the cold season are also found to be 
externally forced, and larger than could be as a result of natural 
(internal) variability alone (P < 0.05). Such changes enforce heat- 
trapping systems, and serve to systematically enhance atmo- 
spheric stability in cold season together with decreases in winter- 
time low-cloud’s cooling effect, supporting the configuration of 
the warming pool. 


The evolution of anthropogenic aerosols plays an important 
role in the northeast Pacific SST trends subsequently (post-1996) 
when East Asian aerosol emissions decreases. The cooling to 
warming SST transition after around 1995/1996 is reproduced in 
the aerosol-forced historical simulation, when the ALL-forcing 
trend pattern is characterized by amplified warming, with the 
contribution of both GHG-induced radiative warming and 
AERindust-induced local warming to the SST increase over the 
region. Our bivariate attribution analysis over 1996-2021, based 
on multiple model data sources, demonstrates that forcing by 
elevated GHG levels and the recent industrial aerosol-load 
decrease are key causes for the here-detected warming pool 
(P< 0.05). These results further strengthen the MHWs attribution 
findings, that GHG forcing has virtually certainly caused the 
2019-2021 MHW, by demonstrating that GHG forcing has also 
discernibly changed the background state against which MHW 
events occur. This indicates that the region will very likely 
experience more intense MHW events as the GHGs emissions 
continue to rise. 

The here-detected long-term Pacific warming pool together 
with the strengthening high-pressure system on the background 
state of the high-impact MHWs, which we show that significantly 
deviate from natural variability, and that will continue and 
intensify in the course of unfolding anthropogenic climate 
change, will provide favorable conditions over the northeast 
Pacific for even more sever marine heatwave events in the future. 
Such change could have profound societal and marine ecosystem 
impacts over the region. 


Methods 

Observations. The area of research is the North Pacific ocean defined by the 
10°N-60°N, 90°E-90°W latitude-longitude domain. Two datasets of sea surface 
temperature are used: (1) NOAA OISSTv2”, which is daily remotely sensed 
National Oceanic and Atmospheric Administration (NOAA) Optimum Inter- 
polation (OI) SST V2 high-resolution (0.25°) gridded SST data available over 
January 1982 to April 2021. (2) The Hadley Centre HadISST®’, which is produced 
monthly with 1° grid resolution since 1870. The mean sea-level pressure (SLP) and 
500 hPa geopotential height (500Gph) are derived from ERA5“? reanalysis. Cloud- 
cover data is from EUMETSAT’s Satellite Application Facility on Climate Mon- 
itoring (CM SAF) available over 1982-20182. The summary of observation and 
model data used in this study is presented in Table 2 and Supplementary Table 1. 


Large ensembles. Two ensembles of ocean-atmosphere coupled model simula- 
tions are used; the Community Earth System Model 35-members Large Ensemble 
(CESM-LE)28, and the Max Planck Institute for Meteorology 100-members Grand 
Ensemble (MPI-GE)°°. The CESM-LE includes 35 simulations (members) running 
from 1920 to 2100. From 2006 to 2100, the Representative Concentration Pathway 
8.5 forcing (RCP8.5%) is used. The MPI-GE is conducted by using the Max Planck 
Institute Earth System Model (MPI-ESM1.1) and includes 100 simulations 
(members) running from 1850 to 2100 under the RCP8.5 scenarios. For disen- 
tangling the forced response and the internal climate variability each member in 
the ensembles is initialized with different initial conditions. The large size of the 
ensemble is a crucial requirement to robustly sample internal variability. The mean 
of each ensemble averages out the internal variability, thus represents the forced 
response of the system. 

In addition, for attributing the detectable signal to a specific forcing agent we 
make use of two ensembles that are identical to the CESM-LE”%, but each ensemble 
excludes the time evolution of one forcing agent: greenhouse gases (LE-fixGHG, 
include 20 members) and anthropogenic aerosols (LE-fixAERindust: include 20 
members). 


CMIP6-DAMIP. We also analyze five models that participate in the Detection and 
Attribution Model Intercomparison Project (DAMIP**) component of the Coupled 
Model Intercomparison Project Phase 6 (CMIP6*°). From the DAMIP** single- 
forcing runs, we consider two groups of simulations. One group (GHG) includes 
34 simulations conducted by five models forced with historical well-mixed 
greenhouse gases only. A second group (AER) includes 34 simulations conducted 
by 5 models forced with anthropogenic aerosols only. The DAMIP models 
(ensemble size) are CNRM-CM6-1 (6), CanESM5 (10), GISS-E2-1-G (4), 
HadGEM3-GC3-LL (4), and IPSL-CM6A-LR (10). In these simulations, aerosol 
and GHG emissions (or concentrations) are allowed to vary in time whereas all 
other forcing variables are set to pre-industrial values. The formula to give equal 
£ ~ where d is the number of models 


intl 


weighting of the individual models is, n = 


COMMUNICATIONS EARTH & ENVIRONMENT | (2022)3:131 | https://doi.org/10.1038/s43247-022-00461-2 | www.nature.com/commsenv 9 


ARTICLE 


COMMUNICATIONS EARTH & ENVIRONMENT | https://doi.org/10.1038/s43247-022-00461-2 


Table 2 Observation and model data used in this study. 


Source 


NOAA OISST42 V2; HadISST&8 

ERA5“? reanalysis (1979-2021) 

ERAS reanalysis (1979-2021) 
EUMETSAT®2 satellite data (1983-2018) 


Observation data 


Sea surface temperature (SST) 
Mean sea-level pressure 
Geopotential height at 500 hPa 
Cloud-cover fraction 


Model data 


Historical simulations in 
monthly timescale (ALL forcing) 
Historical simulations in daily 
timescale 

GHG-forcing only, AER-forcing 
only runs 


Source 
CESM1-LE2? (35) and MPI-GE2° (100) 


CESM1-LE2? (20 members; ALL and fixGHG) 


34 runs conducted by 5 models of 
CMIP635-DAMIP34 

20 runs conducted by CESM1-LE2? 
Pre-industrial control CMIP635 
simulations (7000 years) 
Paleo simulations (0850-1850 


millennium) 


CCSM4°° Paleo simulations 


and | is the ensemble size”. The final internal variance is then just 1/n the internal 
variance. Thus, in the multimodel ensemble mean of 34 simulations conducted by 
five models (GHG-only and AER-only simulations), the internal variability is 
reduced by about 70%, which leads to an enhanced signal-to-noise ratio in esti- 
mated signal patterns. The name of the models from DAMIP-CMIP6, the number 
of model ensemble members and the pre-industrial control years used in the study 
is presented in Table 2 and Supplementary Table 1. 

The summary of single-forcing experiments used in this study is as follows: 


e@ ALL signal: Ensemble of 35 runs from CESM-LE, and the ensemble of 100 
runs from MPI-GE forced with ALL forcing, which includes anthropogenic 
factors such as human emissions of greenhouse gases, atmospheric 
aerosols, ozone, land-use changes, and natural external factors such as 
stratospheric aerosols due to the large volcanic eruptions and solar forcing. 

@ GHG signal: Ensemble of 20 runs from CESM-LE, and the ensemble of 34 
runs conducted by 5 CMIP6-DAMIP models forced with historical changes 
in well-mixed greenhouse gases. 

e AER signal: Ensemble of 20 runs from CESM-LE, and the ensemble of 34 
runs conducted by 5 CMIP6-DAMIP models forced with historical changes 
in anthropogenic aerosols. 


Defining marine heatwaves. We identify marine heatwaves (MHWs) from daily 
NOAA OISST time series available from January 1982 to March 2022 and follow 
the standardized and widely used!-!?,7! MHWs definition developed in Hobday 
et al.!°. MHWs occur when SSTs exceed a seasonally varying threshold, defined as 
the 95th-%tile of SST variations based on a 30-year climatological period 
(1983-2012), for at least five consecutive days. At each location and for each 
MHW, we calculated the event duration (time between start and end dates), fre- 
quency and intensity (SST anomaly above the threshold average over the event 
duration). In the analysis, two events with a peak of less than 3 days are considered 
as a single event. The MHW definition used in this study is available as software 
modules in R (heatwaveR”2). 


Extreme event attribution. Extreme event attribution*? is used to assess the 
influence of GHG forcing on the duration and intensity of the observed MHWs. 
Following Hannart et al.*°-4!, we present event attribution in terms of necessary 
and sufficient causation. In order to obtain reliable estimates of the probabilities, 
we use daily SST output from CESM1-LE with a 20-member ensemble with ALL 
forcing and a 20-member ensemble with excluded time evolution of GHG forcing 
(LE-fixGHG). The differences among 20 ensemble members are due to internal 
variability and the 20 simulations can be considered as 20 plausible realizations of 
the real world?!. To the best of our knowledge, the CESM1-LE is the only com- 
prehensive model available with complementary historical single-forcing large 
ensembles in daily timescale. Given that the magnitude of daily SST variability in 
CESM1-LE compares well with observations (standard deviation of 1.91 vs 1.94 °C, 
respectively, based on detrended data during 1982-2021), this model ensemble can 
be used for MHWs attribution analysis. In addition, the results of detection and 
attribution analysis of long-term SST presented in “Detection of systematic changes 
in the Pacific warming pool” and “Univariate trend detection and bivariate trend 
attribution”, indicate that the CESM-LE is suitable for an analysis of the attribution 
of MHWs in the region because of its strong detection of the anthropogenic signal 
(scaling factor a; is very close to unity; Figs. 4e and 6a). 

Choosing the 1982-2021 time period, the MHW characteristics (duration and 
intensity) from each year and each 20 realizations of CESM1 model are pooled 
together (780 years total) and probabilities are estimated. For each detected MHW, 
we estimate the probabilities that a MHW has occurred that equals or exceeds the 
duration and intensity of the observed MHW in ALL forcing (actual world) and 
fixGHG forcing (counterfactual world) scenarios. 


We calculate the probability of necessary causation (PN), sufficient causation 
(PS) and probability ratio (PR) separately for MHWs duration as: 


duration duration 
pauration pauratior | — paura 
ALL fixGHG , ALL 
PR= piuration ? PN=1 pduration ? PS=1 | — Pduration * Q) 
‘fixGHG ALL ‘fixGHG 
and similarly for MHWs intensity: 
tensity intensity intensity 
pin mst) 1-P 
ALL. JixGHG ALL 
PR= pintensity ; PN=1 intensity ? PS=1 = pintensity (2) 
fixGHG ALL fixGHG 


It should be noted that the equations presented here only apply as PN or PS if 
the resulting values are greater than 0; if negative, the PN or PS is assigned a 
probability of 0. 


Calculating uncertainties on PN. To calculate the uncertainty on the probability of 
necessary causation (PN) a resampling*? method is used. The pool of data from 
each year in 1982-2021 and from each 20 ensemble member in CESM1-LE is re- 
sampled to reproduce a set of data to use for the calculation of the probabilities. 
With 1000 estimates of the density curves, a nonparametric 90% confidence 
interval for each of the metrics can be determined. 


Estimating natural (internal) variability. We used three sources for the estima- 
tion of time-invariant and time-evolving (historical) internal variability of climate 
system, of which PDO, ESNO, NPGO, etc. are part, as well as natural (internal + 
external) variability. 


1: Time-invariant internal variability. The long pre-industrial control simulations 
from global climate models (GCMs) participating in the CMIP6 project are per- 
formed under control conditions (i.e., with constant atmospheric composition, no 
episodic volcanic influences, and no variation in solar output). The models 
(number of years used from control integration) are CESM1-LE (1000), MPI-GE 
(2000), CNRM-CM6-1 (500), CanESM5 (1000), GISS-E2-1-G (1000), HadGEM3- 
GC3-LL (500) and IPSL-CM6A-LR (1000). Since climate models may under- 
estimate variability, we double the simulated variance prior to the attribution 
analysis’>. The 7000-year pre-industrial control (PIC) runs, which are the con- 
catenated PIC runs of 7 models, provide up to 280 pseudorealizations (the series is 
split into 280 nonoverlapping 25-year segments) of how the climate might have 
changed in the absence of external influences (""Pre-industrial variability"). The pre- 
industrial variability are calculated using concatenated control runs in order to 
capture the very low-frequency 10-20 years fluctuation of SST. For all piControl 
simulations, a linear trend is subtracted, to reduce a possible tiny influence of 
model drift. 


2: Time-evolving (historical) internal variability. We further utilize MPI-GE”° large 
initial-condition ensemble, which produces 100 realizations of single model’s 
response to ALL forcing. This ensemble can provide better estimates of an indi- 
vidual model’s historical internal variability and response to external forcing>!*2. 
The transient forced response (F,) is quantified by taking the ensemble mean of 100 


members at each time step, F, = sie! “, where f,; is a single ensemble member at 
time step t. The “evolving internal variability” is estimated at each time step by 
removing the ensemble mean (F;) from each 100 ensemble members (f,;) and 


calculate the standard deviation across the ensemble”. 


3: Natural (internal + external) variability. We use the Paleo simulations over the 
0850-1850 millennium derived from CCSM45° model to obtain an estimate of 
natural external variability of SST, associated with solar variability and strato- 
spheric aerosols due to large volcanic eruptions. This millennium simulation is split 
into 39 nonoverlapping 25-year segments, which provides 39 pseudorealizations of 
how SST might have changed in the absence of anthropogenic influences*°. 

With the estimated internal variability, derived from the three sources 
mentioned above, we test the null hypothesis that the observed trend over 
1996-2021 is within the 2.5-97.5%-tile distribution of unforced trends (as derived 
either from the pre-industrial control simulations or MPI-GE 100-member 
historical runs) or naturally forced trends (as derived from the 850-1850 
millennium simulation). If the null hypothesis is rejected, it indicates that the 
observed SST changes deviates significantly from internal variability, that is, the 
observed change in SST cannot be explained by internal variability alone, and 
externally (systematically) forced changes are detectable in observed SST record at 
5% significant level. We note here the adoption of risk of false rejection (P< 0.05) 
of the null hypothesis of “no external forcing”. Since when the regional null 
hypothesis is valid, on an average n = 0.05 m local alternatives will falsely be 
rejected (m number of grid points)”*. 


Long-term trend detection and attribution. We follow the same detection and 
attribution approach used in previous studies by Barkhordarian et al.373855-67 that is 
a two-step process. In the first step—detection—consists of assessing whether sys- 
tematically (externally) forced changes are detectable in the observed SST change. 
This is achieved by testing the null hypothesis that the observed change in SST is 
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drawn from the population of an undisturbed climate. In the second step—attribution 
—we examine the null hypothesis that the observed change in SST is drawn from a 
hypothetical population of a climate disturbed by a specific external influence. 

The attribution analysis here is based on estimating the amplitude of the 
response of SST to each external forcing from the observations via the estimation of 
scaling factors’*’°, which is a linear regression model as follows: 


Yobs = UC; — Uj); + Ups (3) 


where y,»; represents the observations and each x; the modeled response to one of 
m forcings that is anticipated by climate models. a; is an unknown scaling factor. 
The noise on y,,;, denoted by u,,;, is assumed to represent internal climate 
variability, while the noise on x; denoted by u;, is a result of both internal 
variability and the finite ensemble used to estimate the model response. In order to 
account for the noise in response patterns Total Least Squares (TLS) 
methodology’’, which minimizes the perpendicular distance between the scatter 
points and the best-fit line is used. 

Uncertainties in a; are estimated by accounting for the effect of internal climate 
variability on y,,;, using samples from climate model control simulations. 
Therefore, the distribution of the scaling factor a; is assessed from fits of the 
regression model (Eq. (3)) to 280 nonoverlapping control run segments derived 
from 7000-year control simulations (conducted by 7 global climate models). 

Detection of a climate change signal occurs if the uncertainty range around a 
scaling factor a; is shown to be significantly different from zero. This is handled by 
testing the null hypothesis Hp;: a = 0 (where 0 is a vector of zeros). If the null 
hypothesis Hp, is rejected, it indicates that the observed representation of climate 
change deviates significantly from internal variability, that is, the observed change 
Yovs cannot be explained by internal variability u,,, alone. Once detection has been 
established, attribution (consistency of observed changes with a combination of 
external forcing) is assessed by testing the null hypothesis H47: a = 1 (where 1 
denotes a vector of unit). When there is insufficient evidence to reject H47, the 
attribution of changes to the respective forcing is claimed>7-38.67, 


Data availability 

The NOAA OISST dataset is publicly available at https://downloads.psl.noaa.gov/ 
Datasets/noaa.oisst.v2.highres. The climate model simulations are available via the Earth 
System Grid Federation (ESGF) archive of Coupled Model Intercomparison Project 6 
(CMIP6) data (https://esgf-index1.ceda.ac.uk/projects/esgf-ceda/). The ERAS reanalysis 
data can be obtained from https://cds.climate.copernicus.eu, and the observed cloud- 
cover data from https://www.eumetsat.int. The CESM1 Large Ensemble (LE) data are 
publicly available at https://www.earthsystemgrid.org/dataset. 


Code availability 
The marine heatwave definition used in this study is available as software modules in R 
(heatwaveR7). 
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